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SUMMARY 

A direct transfer matrix method is developed to analyze the 
linear and nonlinear dynamics of multiple-load-path bearingless 
rotor blades. The method is applied to determine (1) the natural 
frequencies and modes about the initial state, (2) the nonlinear 
steady state deflections, (3) the natural frequencies and modes about 
the steady deformed state and (4) aeroelastic stability of multipie- 
load-path rotor blades in hover. A Newton-Raphson iterative 
method based on quasilinearization of the nonlinear distributed 
boundary value problem is developed to solve the steady state 
deflections of the blade. Aerodynamic forces are calculated 
employing two dimensional strip theory and quasi-steady 
aerodynamics. The formulation is validated by comparing the results 
for single and multiple-load-path blades with those obtained by 
other methods in the literature. For forward flight applications a 
discretization based on either modal coordinates or harmonic 
analysis is recommended. 
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1. INTRODUCTION 

A large class of systems occurring in engineering practice 
consist of one dimensional beam and beamlike elements. Sometimes 
the systems contain either a single element or a number of elements 
linked together end to end in the form of a chain. Well known 
examples are simple beams, rotor or propeller blades, continuous 
beams, fuselage bulkheads, and turbine generator shafts. The 
transfer matrix (which is one form of frequency response matrix) is 
ideally suited to treat such one dimensional chainlike systems 
governed by linear equations. Intermediate conditions and the 
number of degrees of freedom have no effect on the order of the 
transfer matrices and depends only on the order of the governing 
differential equations. 

Holzer [1] initially applied the transfer matrix method to 
determine the torsional vibrations of rods and the method is 
generally known as Holzer's method. Myklestad [2] applied a method 
analogous to Holzer's method to determine the bending-torsion 
modes of beams and the method is usually called Myklestad's 
method. Thomson [3] applied the method in a matrix form to more 
general vibration problems. The original application of the transfer 
matrix method also includes the description of steady state 
behaviour of four terminal electrical networks, in which case the 
method is commonly referred to as four— pole parameters method. 
Molloy [4] applied four pole parameters to study acoustical, 
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mechanical, and electromechanical vibration problems. Pestel and 
Leckie [5] have listed transfer matrices for elasto-mechanical 
elements up to twelfth order and the textbook contains several 
references on transfer matrices. Rubin [6,7] has provided a general 
treatment for transfer matrices and their relation to other forms of 
frequency response matrices. Transfer matrices have been applied 
to a wide variety of engineering problems by a number of 
researchers, including Targoff [8], Lin [9], Lin and McDaniel [10], 
Mercer and Seavey [11] Mead [12], Mead and Sen Gupta [13], 
Henderson and McDaniel [14], McDaniel [15], McDaniel and Logan 
[16], Murthy and Nigam [17] a|^d Murthy and McDaniel [18,19]. 
These applications deal with beams, beam type periodic structures, 
cylindrical shells and stiffened rings. 

The main advantage of the transfer matrix method is the 
smallness of the order of the matrices involved. The order of the 
transfer matrix will be equal to the number of elements in the state 
vector. The simplicity gives rise to several numerical difficulties in 
using transfer matrices. These can occur first, when calculating 
higher natural frequencies and their associated mode shapes and 
second, when intermediate geometric compatibility conditions are 
stiff. Despite these numerical problems transfer matrices offer an 
efficient means to study the dynamics of one dimensional systems. 
Combined with leaps in computing power, transfer matrices make it 
possible to tackle new classes of problems such as near real-time 
simulation and optimization. These new problems are feasible 
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because one dimensional systems can be modeled with 
computationally efficient transfer matrices. 

The transfer matrix method is very popular for the analysis of 
rotor blades [20-26] and the reasons are (1) most of the time, only a 
few number of lower natural frequencies and their associated mode 
shapes are of interest, (2) no intermediate stiff conditions are 
involved in the calculation of frequencies, (3) the order of the 
frequency determinant is at the most six by six and (4) the method is 
very appealing for programming. In fact, several rotor dynamics 
programs within the helicopter organizations employ transfer matrix 
modeling for their blades. Examples of such programs include 
Myklestad program of Bell Helicopter [27] , Rotorcraft Airframe 
Comprehensive Aeroelastic Program (RACAP) of McDonnel Douglas 
Helicopter Co. [28], G400 program of United Technologies [29] , C60 
Program of Boeing Helicopter and KTRAN of Sikorsky Aircraft. 

All the existing helicopters employ either articulated, teetering, 
gimbaled or hingeless blades for their rotors to relieve high blade 
stresses encountered at the blade root during normal operating 
conditions. Examples of these rotor types are shown in Figures 1.1 
through 1.4. An advanced configuration known as the bearingless 
rotor is currently being employed in new helicopters. The 
bearingless design is an attempt to realize the several best features 
of articulated rotors (lower vibration and gust sensitivity), teetering 
rotors (low cost) and hingeless rotors (high control power, superior 
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Figure 1.1 Model of An Articulated Rotor 
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Figure 1.3 Gimballed Rotor Hub 




Figure 1.4 Model of a Hingeless Rotor 








8 


flying qualities, mechanical simplicity and low maintenance). The 
recent advances in composite materials make the bearingless rotor a 
viable concept [30]. The rotor blades in these systems are attached 
to the hub via composite flexible structural elements known as load 
paths or branches. These branches accommodate rotor motions by 
deflecting elastically thereby eliminating the need for hinges. 

Conventional helicopter blades are idealized as single-load-path 
blades and as mentioned earlier the transfer matrix method is a very 
convenient and efficient method to analyze such single branch 
structures. All practical designs of bearingless rotors include 
multiple branches at the root, and the one that was flight tested by 
Boeing Helicopter Co. has three load branches consisting of two 
flexbeams and a nonenclosing torque tube [31]. It is shown in Figure 
1.5. Now, almost all of the current bearingless rotor designs consist 
of a single flexbeam with a wrap-around torque tube called a pitch 
cuff (Figure 1.6). 

The goal of this research is to develop the transfer matrix 
method to treat nonlinear autonomous boundary value problems 
with multiple branches. The application is the complete nonlinear 
aeroelastic analysis of multiple-branched rotor blades. Once the 
development is complete, it can be incorporated into the existing 
transfer matrix analyses mentioned previously. There are several 
difficulties to be overcome in reaching this objective. The 
conventional transfer matrix method is limited in that it is applicable 
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Figure 1.5 Components of a First Generation Bearingless Rotor 
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Typical Current Bearingless Rotor Design 
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only to linear single branch chain-like structures, but consideration 
of multiple branch modeling is important for bearingless rotors [33]. 
Also, hingeless and bearingless rotor blade dynamic characteristics 
(particularly their aeroelasticity problems) are inherently nonlinear. 
Murthy and Joshi [34] have developed a transfer matrix method to 
determine the natural frequencies and modes of rotor blades with 
multiple branches at the root. This development was based on linear 
homogeneous equations of motion. Sangha [35] has described a 
transfer matrix method to analyze a bearingless multiple-load-path 
blade. 

In the present work the nonlinear equations of motion and the 
multiple-branched boundary value problem are treated together 
using a direct transfer matrix method. First, the formulation is 
applied to a nonlinear single-branch blade to validate the nonlinear 
portion of the formulation. The nonlinear system of equations is 
iteratively solved using a form of Newton-Raphson iteration scheme 
developed for differential equations of continuous systems. The 
formulation is then applied to determine the nonlinear steady state 
trim and aeroelastic stability of a rotor blade in hover with two 
branches at the root. A comprehensive computer program is 
developed, and is used to obtain numerical results for the (1) free 
vibration, (2) nonlinearly deformed steady state, (3) free vibration 
about the nonlinearly deformed steady state and (4) aeroelastic 
stability tasks. The numerical results obtained by the present 
method agree with results from other methods. 
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2 . EQUATIONS OF MOTION 

2.1 Nonlinear Equations for Elastic Bending and Inertial 

Loadings 

The nonlinear differential equations of motion for the fully 
coupled elastic flapwise bending, chordwise bending, torsion and 
axial extension of twisted nonuniform rotor blades are given below. 
The development of the equations is the subject of reference [36], in 
which their complete derivation is presented. In addition to the 
present project, these equations form the basis for structural beam 
modeling in the current state of the art rotorcraft analysis programs. 
The equations are valid for any beam, and the mass, elastic and 
tension axes need not be coincident. The coordinate system they are 
derived in is an undeformed coordinate system, with x positive 
outward along the span, y positive towards the leading edge, and z 
positive upwards (see Figure 2.1). For algebraic conciseness the 
terms containing eA, B*, B 2 , Ci and C* are treated as zero. It should 
be noted that this assumption does not affect the general nature of 
the formulation presented here-in. 


Axial Extension: 

[EA(u'+v' 2 /2+w ,2 /2)]' + Q 2 mu - mu + 2£2mv = -Q 2 mx (2.1) 
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Lead-Lag: 

- (Tv')' + {(EI z cos2(e+<(>) + EIysin 2 (0+<J>)]v" 
+(EI z -EIy)cos(0+<|>)sin(0+<|>)w" } " 

+ 2Qmii + mv - me<j>sin0 - 2mei2(v'cos0+w'sin0) 

- mQ 2 [v+ecos(0+<f>)] - 2mi2pp C w 

- {me[£2 2 xcos(0+<|)) + 2Qvcos0)}' = Ly 


Flap: 

- (Tw')' + {(EI z -EI y )cos(0+<|))sin(0+()))]v" 
+(EI z sin 2 (0+<|)) + EI y cos 2 (0+<|>)w"} 
+ 2m£2p pc v + mw + me<j)cos0 
- {me[Q 2 xsin(0+<})) 

+ 2Qvsin0)}' = L w - Q 2 mP pc x 


Torsion: 

- (GJ+TkJ)*']' + (EI z -EI y )(w" 2 -v" 2 )cos0sin0 

+ v"w"cos20] + + mQ^k^-k^ i )cos20 

- me[fl 2 x(w'cos0-v’sin0) 

- (v -Q 2 v)sin0+ wcos0] = Me 

- m0 2 (k^ 2 -k^ j)cos0sin0 - meQ 2 P pc xcos0 


( 2 . 2 ) 


(2.3) 


(2.4) 
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2.2 First Order Equations 

Equations (2.1)-(2.4) can be reduced to a system of first-order 
nonlinear differential equations in terms of the state vector contain- 
ing deflections and forces [37]. The state vector {Z} is defined as 

L u w v e £ <() M x My M z V y V z V x J (2.5) 


and the governing equations are written in the undeformed 
coordinate system. 

First the differential equations of the three translations u, w 
and v are developed. From equation (90) of reference [36] 


V x = EA 


u 



( 2 . 6 ) 


The bending rotations (slopes) are 

w' = e (2.7) 

v' = $ (2.8) 

Substitution of equations (2.7) and (2.8) into equation (2.6) gives 

u' = -e 2 /2 - 4,2/2 + V X /EA (2.9) 

Expressions for the bending rotations (e\ 4', <!>’) are determined 
as follows. From equations (91), (92) and (96) of reference [36] 


M x ' = GJ<t)' 


( 2 . 10 ) 
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M z ’ = EI z [v" cos(0+())) + w" sin(0-+-<t>)] 
My’ = EIy[v" sin(0+<|)) - w" cos(0+0)] 

In matrix form, the above equations are 


( 2 . 11 ) 

( 2 . 12 ) 


' M x ’ ' 

= 

1 

o 

o 

Q 

1 


' <J>’ ' 

M z ’ 

> = 

EI z sa EI z ca 0 

< 

S’ 

k My 1 > 

= 

_ -Elyca Elysa 0 _ 


. 0’ > 


(2.13) 


where s = sine, c = cosine and a = 0 + <J>. In transforming these 
moments to the undeformed coordinate system, the following 
parameters must be considered. From equation (A6) of [36] 


0 + $ + v'w' = 0 + <j> - J v"w'dx + v'w' 


o 


or 0 + <J> + v'w' = 0 + <J) + v 
where 

x x 

v = v'w' - J v"w'dx = J v'w"dx 
o o 

Assuming small o, 

0 + $ + v'w’ = 0 + <|) 


(2.14) 


(2.15) 


(2.16) 


Given equation (2.16), equations (3) of [36] relate the bend- 
ing moments in the deformed coordinate system (x\ y', z') to those 
in the undeformed coordinate system (x,y,z) 
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M x = M x ' - M z ’(-£sa+ecoc) - My’(^ca+esa) 

(2.17) 

M z = M x ’e + M z 'ca + M y ’sa 

(2.18) 

M y = M x '^ - M z 'sa + My’ca 

(2.19) 


by neglecting e 2 , ^2 compared to unity. In matrix form the above 
equations are 


l 


M x 

My 

M z 


1 

+esa - eca 

^ca - esa 

e 

ca 

sa 

% 

-sa 

ca 



( 2 . 20 ) 


Premultiplication of equation (2.13) by the coefficient matrix of 
(2.20) results in an expression relating derivatives of rotations to 
bending moments in the deformed coordinate system. Substituting 
equation (2.20) into this expression gives 


M x ’ 


' € ' 

M z 

> = ter 

%' > 

l My > 


, <t>’ j 


where [C] is defined as 


b 2 + bi£ -bi^ - b 3 e GJ 


[C] = 


and 


-bi 

-b2 


b 3 

bi 


GJe 
GJ£ _ 
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bi = (EI y -EI z )casa 
b2 = EI y c 2 a + EI z s 2 a 
b 3 = EI y s 2 a + EI z c 2 a 
b4 = b3 - b2 


( 2 . 22 ) 


Inversion of the coefficient matrix [C] in equation (21) and subse- 
quent rearrangement yield the desired differential equations of 
rotations. 


v| /' = GJ(b 3 ^-bie)M x /D + GJbi(l+S 2 +b 3 Se/bi)M z /D 

- GJb 3 (l+e 2 +bi^e/b 3 )M y /D (2.23) 


% = GJ(bi^-b 2 e)M x /D + GJb 2 (l+^ 2 +b^e/b 2 )M z /D 

- GJbi(l+e2+b2^e/bi)M y /D (2.24) 

(J)’ = (b3b2-b 2 1 )M x /D + b3b2-t? 1 )eM z /D 

+ (b 3 b 2 -l? 1 )^M x /D (2.25) 

where 

D = GJ[(b 3 b 2 -l? 1 )(l+4 2 +e 2 )] 

and 

b 3 b 2 -b 2 1 - EI y EI z 


In equations (2.23) through (2.25), are neglected compared to 

unity, yielding 


b3^-bie 


Mv + 


bi 


M 7 - 


b3 w 

ei v ei 7 M y 


EI y EI z 


EIyEI z 


(2.26) 
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<t>’ = 


bi^-b 2 e 

EI y EI z 


M x + 


Mx Mz ^ 
GJ GJ e 


b 2 

EIyEIz 


Mz - 



% 


bi 

EIyEIz 


Mz 


Assuming the torsional displacement angle <J> is small, 
cos(0 + <)>) = cos0 - <j>sin0 

► 

sin(0+<()) = sin0 + <{>cos0 


(2.27) 

(2.28) 


(2.29) 


Substitution of (2.29) into (2.22) gives 


bi = (EI y -EI z )(cos0 sin0 + <j>cos20) 
b 2 = EI y (cos 2 0 - <|)sin20) + Elz(sin 2 0 + <|>sin20) 
b 3 = EI y (sin 2 0 + <(>sin20) + EI z (cos 2 0 + <J>si n20) 
b4 = b3 - b 2 


> 


j 


(2.30) 


neglecting <J> 2 terms. 
Define 


ai = (EIy-EI z )cos0sin9/EIyEI z " 
a 2 = cos 2 0/EI z + sin 2 0/EI y ^ 
a 3 = sin 2 0/EI z + cos 2 0/EI y 
a 4 = a 3 - a 2 ' 

Then from equation (2.30) and (2.31) let 


(2.31) 



20 


b i/EIyEI z = ai - a4<|> 
b2/EI y EI z = a2 - 2a i <(> 
b 3 /EI y EI z = a3 + 2a i 
b4/EI y EI z = a4 

Substituting equation (2.32) into the relations (2.26), (2.27) and 
(2.28) for e', t,' and <>', and neglecting third order nonlinear terms 
yields 

e’ = (a 3 ^-a 1 e)M x + (a r a 4 <|>)M z - (a 3 +2a 1 <())M y (2.33) 

S’ = (a^-a 2 e)M x + (a 2 -2a 1 <J>)M z - (a r a 4 <|>)M y (2.34) 

f = M X /GJ + M z e/GJ + My^GJ (2.35) 

The first order differential equations for the moments are 
developed from equations of equilibrium (68) of inertial moments 
(88) of [36]. 

dM x 

( j x = V y e - V z £ - m{e[(v -Q 2 v)(sin0+<()cos0) 

- w(cos0+<J>sin0) + 2Qusin0 - Q 2 xp pc cos0] . j.2^ 

■ f2 2 (k^ 2 -k^)(cos0sin0) + <|)cos20 - <|> 2 cos0sin0) 

- 2 ^[( k m 2 ’ k m ! )^sin0cos0 

+ e( k2 2 sin20 + k^cos 2 ©)]} - M e (2.36) 

where M e is the aerodynamic pitching moment about the elastic axis 



dM z 

dx = -V y + V x l; + me[Q 2 x(cos0-(f)sin0) + 2£2vcos0] 


(2.37) 



= V z + V x e - me[Q2x(sin0+<|)cos0) + 2Qvsin0] 
ci x 
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(2.38) 


The first order differential equations for the forces are 
developed from equations of equilibrium (69) of [36] and inertial 
forces given by equation (87) of [36]. 

= m[v - e4>sin0 - ft 2 [v + e(cos0-<]>sin0)] 

+ 2Q[u - e{£cos0+esin0)]} - L v - 2mf2P pc w (2.39) 

where Ly is the aerodynamic force acting in the y direction 
positive towards the leading edge 

= m(w - e$cos0) - L w + 2mQp pc v + mP pc Q 2 x (2.40) 

where L w is the aerodynamic force acting in the z (upward) 
direction. 

dVv 

= m0 2 (x+u) - 2ilmv + mii (2.41) 

Thus equations (2.9), (2.7), (2.8) and (2.33) through (2.41) are the 
desired first order nonlinear differential equations governing the 
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state vector. The aerodynamic loadings, L v , L w and M e are evaluated 
next. 


2.3 Aerodynamics 

Helicopter rotors typically operate in a complex aerodynamic 
environment. During one revolution (typically 1/5 of a second) the 
blade can experience large variations in angle of attack and Mach 
number, unsteady aerodynamic effects and stall effects. 

To model the rotor aerodynamic environment, lifting line or 
lifting surface theories can be employed. Lifting line technology is a 
more mature technology and for that reason is used extensively in 
the helicopter industry. 

The rotating blade is modeled as a rotating lifting line. Strip 
theory is then applied to this lifting line, resulting in a finite number 
of discrete bound vortex segments representing airfoil sections. 
Unsteady aerodynamic theory has been developed for the lifting line 
approximation of a fixed wing in references [37] and [39]. Unsteady 
aerodynamics takes the complete vortex system into account and 
gives rise to two classes of forces - circulatory and noncirculatory. 
The circulatory forces are a consequence of the vorticity (bound and 
trailed). The noncirculatory forces (sometimes referred to as virtual 
or apparent forces) are not due to vorticity. Lift may be expressed 
as the sum of its circulatory and noncirculatory components. 


LTotal - ^circulatory + Lfsjoncirculatory 


(2.42) 
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Theodorsen noted in [37] that if a wing was undergoing pure 
harmonic motion at frequency co, the wake vorticity was periodic and 
subsequently circulatory lift would be also. The Theodorsen lift 
deficiency function C(k) thus modifies the circulatory lift. It accounts 
for the effect of shed vorticity in the wake and subsequent time 
varying lift build up on the wing. 

LTotal = C(k)Lcirculatory + LNoncirculatory (2.43) 

where k is a reduced frequency, given by 
cob 

k = — (2.44) 

The magnitude of C ranges from .5 at high frequency to 1. at low 
frequency, and thus usually reduces the circulatory lift. The phase of 
C represents the delay (or lag) in the lift build up. The magnitude 
and phase of C(k) are shown in Figure 2.2. There are different levels 
of approximations for unsteady aerodynamics, and usually they are 
divided into the following three categories. 

(i) Complete Unsteady Aerodynamics 

The complete theory is employed to represent the unsteady 
effects of the shed vorticity, and includes circulatory and 
noncirculatory forces. 


LTotal — C(k)Lcirculatory + LNoncirculatory 


(2.45) 
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The results obtained correspond to characteristic frequencies of 
approximately 40 cycles per second. 

(ii) Quasi-Unsteady Aerodynamics 

In this instance the noncirculatory terms are neglected. 

LTotal = C(k)Lcirculatory (2.46) 

The theory yields satisfactory results for frequencies between 
5 and 15 cycles/sec, 

(iii) Quasi-Steady Aerodynamics 

I^Total = ^Circulatory (2.47) 

This approximation is generally used at characteristic 
frequencies which are on the order of 1 cycle per second. It 
neglects the effect of wake vortices on the flow field, but 
unsteady rotor motions still produce an unsteady downwash. 

This unsteady downwash is employed to establish the bound 
circulation by satisfying the tangency boundary condition. 

In the present formulation the quasi-steady aerodynamics are 
employed because the aeroelastic instabilities under investigation 
are quite low. (In hover the reduced frequency is small, so C(k) = 1 
and negligible phase lag is a good approximation. As forward speed 
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is increased, the reduced frequency also increases. For increasing k, 
C(k) is a value <1.0 (as can be seen in Figure 2.2) and the accuracy of 
the approximation is somewhat reduced although it is still useful.) 


2.3.1 Airloads 

The quasi-steady lift and moment coefficients about the mid 
chord point are given by [39] as 

C T = 2 I V TTp (2 - 47> 

C c/2 = C + 2 /VS^* ** (2 ' 48) 

where the downwash velocity is 


w(y,t) 


P z a ^ P z a 
pt py 


and 

c q t s = iqs/ P u2b 


c m = mq s /pu2b 2 


(2.49) 


The equation for the meanline of a thin airfoil which is free to 
translate vertically and pitch about the mid-chord is given by (see 
Figure 2.3) 




Figure 2.3. Airfoil in Vertical Translation and Pitch About the 

Mid-Chord 
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z a(y>0 = h c /2 + tan ay (2.50) 

From equation (2.49), assuming small a the downwash velocity is 
given by 

w(y,t) = hc /2 - ay - Ua 

Nondimensionalizing with respect to semi-chord b yields 

w ( y,t) = h c /2 - Ua + bay (2.51) 

Equation (2.51) is then substituted into the expressions for quasi- 
steady lift and moment (2.47) and (2.48). The resulting equations 
may then be integrated with respect to y with the help of the 
integrals in Appendix A. 

no 2ic H • 

c 1 = "u (-11C/2+ Ua + 2«) (2.52) 

C "> q c/2 = U ( -^ 2+Ua) (2 ' 53) 

Traditionally the coefficients are referred to the quarter chord 
location (also the approximate aerodynamic center). The following 
equations relate the airloads of mid-chord and quarter-chord and are 
obtained from the geometry of Figure 2.4. 






1 = 1 c /4 = lc/2 ; a = a C /4 = «c/2 
h c/2 = h c /4 - ab/2 
m c/2 = rn c /4 + lb/2 


30 




(2.54) 


Substitution of equations (2.52) and (2.53) into equation (2.54) and 
noting the thin airfoil theory theoretical result that 'a' (lift curve 
slope) equals 2k yields the two dimensional section airloads 


1 = pUa(-h c / e + Ua + ba 


(2.55) 


b 

m c / 4 = -pbUa (-) 2 a (2.56) 

Consider the airfoil to have properties such that the shear center 
(elastic axis) is not coincident with the quarter-chord location. 
When the airfoil is plunging and pitching about the shear center as 
shown in Figure 2.5, the following relations may be written 


h c /4 = h e + ae b 
m c / 4 = m c - leb 


(2.57) 


Obtain 1 and m e with respect to the elastic axis (or shear center) 
(consistent with the nonlinear equations of motion for elastic 
bending). Substitution of equations (2.55) and (2.56) into equation 
(2.57) and replacement of b by c/2 yields 
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1 = U{-h e + Ua + |- (l-e)a) (2.58) 

me = 9 ~f U(f) 2 (f (-he + Ua) + a(4e-e2- 1 } (2.59) 

Ultimately it is desired to express 1 and m c in terms of elastic 
deflections, and induced and angular velocities. Towards this end, 
consider a two dimensional airfoil in unsteady motion as shown in 
Figure 2.6, noting that U is shown opposite to the free stream 
direction to show it as the airfoil motion. Assuming cosa = 1 and 
sina = a the following expressions for the velocities along the y', z 
axes can be written 

U~Ut (2.60) 

Up = h c - Ua (2.61) 

From Appendix A of reference [40], the velocity components can be 
written as 

Ut = £ix + v (2.62) 

Up = -Qx(0+<|)+'u) - (0+<J>)v + vj + w + £2vx(e+Pp C ) (2.63) 

a + Q(e+P pc ) (2.64) 


where 
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X 

■o=J ^e'dx 

o 


nonlinear angle-of-attack term 
due to bending 


(2.65) 


Knowing the aerodynamic loadings at the correct location along 
the chord line (i.e. at the shear center), resolve the airloads from the 
wind axes back to the undeformed structural axis coordinate system. 
Consider the aerodynamic loadings acting on the airfoil shown in 
Figure 2.7. Resolving loadings into components T and S along the 
deformed (x',y',z') coordinate system axes local to the airfoil gives 


T = lcosa + dsina (2.66) 

S = lsina + dcosa (2.67) 


where d is the profile drag given by 
d = J pU 2 cd 


( 2 . 68 ) 


and 1 is given by equation (2.58). From the geometry of Figure 2.6 


cosa 


= U T = V u p+ u t “U T /U 


(2.69) 


cosa = -Up = ^Up + = Ut/U 


(2.70) 


since 




Figure 2.7. Aerodynamic Forces Resolution 
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_ 

U = a/u£ + Up 


(2.71) 

— 

Substituting the equations for section lift and drag (2.58), (2.68) and 
(2.69) through (2.71) into equations (2.66), (2.67) and (2.59) gives 
airloads in the deformed coordinate system (x\y’,z') as follows. 

— 

pac c 

T = 2 {-UpUt + 2 (l-e)U T a) 


(2.72) 


S = P 2 { -Up - 2 (l-e)Upa ■ -f u£} 


(2.73) 


Assuming that the resultant airfoil velocity V is approximately equal 
to the freestream velocity U in Figure 2.7, a » P + <J> is small. The 
aeodynamic forces may be referred to the undeformed coordinate 
system (x,y,z) used to derive the nonlinear elastic bending equations 
so that they can be incorporated as forcing functions into the 
differential equations of motion of the system. 

L v = S - T((3+<t>) (2.74) 

L w = S(P+<t>) + T (2.75) 

Substitution of equations (2.62) through (2.65), (2.58), (2.72) and 
(2.73) into equations (2.59), (2.74), and (2.75) yields the following 
equations for the aerodynamic forces. 
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pac 9 „ „ 

L v = 2 f v i " £2 2 x 2 k2-£2xvi(0+<|>) 

- {2Qxk2 + (0+<J>)vj } v + {2vj - i2x(0+<J))w 

- k 1 4 > v i — k ii2vj(e+ Ppc) - Qkj w(e+ Ppc)] (2.76) 

L w = [-Qxvi + fi 2 x 2 (0+<|)+'u) - C2 2 xv(e+ P pc ) + Q 2 xki(e+ p pc ) 

- {Qki(e+ p pc ) + £2x(0+<J>) - Vj)v 

- Qx w + k]Qx<|>] (2.77) 

pac c 2 

M e = k 3 L w - ~ 2 ~ {^x<j> + Q 2 x(e+ p pc ) + £2v(e+ p pc ) } (2.78) 


where 

ki = c(l-e)/2 
k2 = Cd/a 
k 3 = ce/2 

The following assumptions are made in the simplification process 
leading to equations (2.76) through (2.78). 

1. All third order nonlinear terms are dropped. 

2. u, v, w, v, <)), u, v, w, e, 0, p pac , vi and k2 are 
assumed to be small. Observe 0 is not assumed to be 
small for elastic and inertial forces, but is assumed to be 
small for the aerodynamics. 

3. All nonlinear terms in u, v, w, e, 4,<|> are dropped. 

The induced velocity is computed as shown below. 
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2.3.2 Induced Velocity Model 

The induced velocity is computed from the combined 
momentum and blade-element theory as given by the following 
equation (reference 41). 


Vi = 


-aa r 

T6” + 



ao r 2 r 
1 6 R 


(e+<j>) 


(2.79) 


Note that it is apparent from equations (2.76) through (2.79) that the 
shear center need not coincide with the aerodynamic center at the 
quarter chord point. 
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3 . BRANCHED BLADES 

Bearingless rotors feature blades connected to a hub via 
multiple branches. Each branch is designed to react a specific 
component of load developed by the rotor blade. Flexbeams are soft 
in torsion, but stiff in bending. They react in and out-of-plane 
bending moments and centrifugal loadings. Torque tubes or pitch 
cuffs (which are aerodynamic fairings enclosing the flexbeams) are 
soft in bending but stiff torsionally; they are used to control the 
blade in pitch. 

A schematic diagram of a multiply branched blade is shown in 
Figure 3.1. The three branches are connected to the blade through a 
rigid clevis. The inboard end of the branches are shown clamped. 
However, the formulation is flexible enough to handle any other root 
and intermediate clevis conditions. 

For the analysis of branched structures by the transfer matrix 
method, it is important to establish the equilibrium and compatibility 
relations across the clevis. 

3.1 Equilibrium Across the Clevis 

The plane of the clevis is shown in Figure 3.2 where h* and h* 

are the y and z locations of the i th branch from the location of the 
blade (point 0). The free-body diagram for this clevis is shown in 
Figure 3.3. Force and moment equilibrium across the clevis yield the 
following equations: 





z 



Figure 3.2. Geometry in the Plane of the Clevis 



2 
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( 3 . 1 ) 


where n = number of branches. 

The above equations can be arranged into a matrix form as 
shown below: 

(f2) = £[Ai] (f>) (3.2) 

i=l 


where 


( f2 } = L M X M Z My Vy V z V X J 
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1 

K 


0 


[I] 1 

0 

0 

-h yi 

[B‘] = 

1 

- - - | 

0 

0 



[0] | 


[I] 



[I] = Identity matrix 
[0] = Null matrix 


3.2 Compatibility Across the Clevis 

Since the clevis is assumed to be rigid, the displacements and 
rotations of the branches and the blade are related as shown below 
(see Figures 3.1 and 3.2) 


u 2 = uj + hje 2 + h ^ 2 

w 2 = w 2 h y '4 


V2 = vj + 

£2 = ej 

4*2 = <t >2 " © r c 




J 


(3.3) 


where 6rc = root collective pitch. The above compatibility equations 
can be arranged into a matrix form as shown below: 


{d 2 } = [B*j {d^} + {b} 


(3.4) 



45 


where 

{d2} T =Lu w v e£<|>J 

{b}T = L o o o o o e rc J 
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4 . NONLINEAR STEADY STATE DEFLECTIONS 

4.1 Equations of Motion 

Substitution of the aerodynamic forcing functions (equations 
(2.76)-(2.79) into the twelve governing first order equations (2.7), 
(2.8), (2.9) and (2.33)-(2.41) yields a set of time dependent nonlinear 
differential equations. The time dependent terms are dropped to 
compute the nonlinear steady state deflections, and the resulting 
equations are given below in matrix form. 

(z(x)}' = [A(x,z)] { z(x) } + {f} (4.1) 

where the state vector consists of deflections, slopes, moments and 
shears. 

{z(x)} T = Lu v w e £ <(> M x M z M y V y V z V x oj 

The nonzero elements of the coefficient matrix are 

A(l,4) = -e/2; A(l,5) = -%H\ A(l,12) = -1/EA; 

A(2,4) = 1; A(3,5) = 1; A(4,4) = -aiM x ; 

A(4,5) = a 3 M x ; A(4,6) = (-a 4 M z - 2aiM y ); A(4,8) = ai; 

A(4,9) = -a 3 ; A(5,4) = -a 2 M x ; A(5,5) = aiM x ; 

A(5,6) = (-2aiM z + a 4 M y ); A(5,8) = a 2 ; A(5,9) = -ai; 

A(6,4) = M Z /GJ; A(6,5) = M y /GJ; A(6,7) = 1/GJ; 

A(7,3) = mQ 2 esin0; A(7,4) = (V y + k 3 k 4 ft 2 xv 
- k 3 kik 4 fl 2 x - k4c 2 H 2 x/16); 
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A(7,5) - -V z ; A(7,6) = m£2 2 (km2-km2)( cos2 0 ' <|>cos0sin0) 

+ Q 2 ecos0v - k 3 k 4 fi 2 x 2 ; 

A(7,13) = k 3 k 4 Q 2 x 2 ; A(8,5) = V x ; A(8,6) = -ma 2 exsin0; 

A(8,10) = -1; A(9,4) = -V x ; A(9,6) = -m£2 2 excos0; 

A(9,ll) = 1; A(10,3) = -mQ 2 ; A(10,4) = k 4 Qviki; 

A(10,6) = mQ 2 esin0 + k 4 Qxvj; A(ll,4) = -k 4 (ki-v)£ 2 2 x; 

A(ll,6) = -k4Q 2 x 2 ; A(ll,13) = -k 4 Q 2 x 2 ; A(13,5) = (M z ai-M y a 3 ) 

where kj = c(l-e)/2 k 3 = ce/2 

^2 = cd/a k 4 = pac/2 

The nonzero elements of the nonhomogeneous excitation vector are 


f(7) - 0 2 m(k^ 2 -km ! )cos0sin0 + k 3 k 4 Qxvi - k 3 k 4 Q 2 x 2 0 

c 2 

k 3 k4Q 2 xkiPp C + k4 j ^ Ppc 

f(8) = Q 2 excos0 
f(9) = -mQ 2 exsin0 

f(10) = -mQ 2 ecos0 + k 4 {-vi + Q 2 x 2 k 2 + nxvj0 + kiQvip pc } 
f(ll) = k 4 Qxvi - k4fl 2 x 2 0 + mQ 2 xp pc - k 4 Q 2 xkip pc 

f(12) = -mfl 2 x 

4.2 Iteration Scheme for Nonlinear Differential Equations 

The overall scheme is to first solve a standard, linear problem 
and then use these results to simplify and solve the more 
complicated nonlinear problem. 
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The following two equations are linearized forms of the 
nonlinear equation (4.1). Both satisfy all the appropriate boundary 
conditions. 

{zO^x)}’ = [A(x,z) z=0 ] {zO)(x}} + {f} (4.1a) 

(z(2)(x)}' = [A(x,z) z=z(1) ] (z(2)(x)) + {f} (4.1b) 

The two solutions denoted by {z(!)(x)J and {z(2)(x)j are defined as the 
first and second approximate solutions of equation (4.1). 

Equation (4.1a) is simply equation (4.1) linearized by the 
starting assumption (z(x)} = {0} in the coefficient matrix [A]. Its 
solution (zO)(x)} = is a linear solution to equation (4.1). 

The solution {zO)(x)} is then used to linearize the coefficient 
matrix [A(x,z) z=z (i)] in equation (4.1b). Solution of equation (4.1b) 

yields (z(2)(x)}, and (z(2)(x)} represents the first nonlinear solution to 
equation (4.1). The second nonlinear solution is based on 

(z(3)(x)} = (z(2)(x)} + (Az(2)(x)} (4.2) 

where (Az(2)(x)} is an incremental solution. To a first order, the 
incremental solution is governed by the following equation. 

{Az(2)(x)}’ = ([A(x,z(2))] + [B(x,z(2))]) {Az(2)(x)} 

- (y ( u 2) (x)} (4.3) 

where (y^(x)} is the unbalanced load and is given by 

{/J} = {z (2) 00}' - [A(x,z( 2 ))] (z( 2 )(x)} - {f} 
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Equation (4.1b) can be used to rewrite the above equation for the 
unbalanced load as 

{y ( u 2) 00} = ([A(x,zO))] - [A(x,z(2))]> (z< 2 >(x)} (4.4) 


The j th column of matrix [B(x,z( 2 ))] denoted by (Bj) is given by 


{Bj} =|V~[A(x,z)] 

I® 2 ! 


=z(2) 


{z(2)} 


/Z=Z 


(4.5) 


where zj is the j th element of the state vector {z}. Since (z( 3 )(x)} and 
{ z(2)(x)} satisfy the same set of boundary conditions, the incremental 
solution vector (Az(2)(x)} satisfies the following zero boundary 
conditions by virtue of equation (4.2). 

(Az(2)(x)} = {0} at boundary points (4.6) 


Now the incremental solution is obtained by solving equations (4.3) 
through (4.6). Once (Az(2)(x)} is known the third approximate 
solution (z(3)(x)} is given by equation (4.2). In general, for the i th 
iterate equations (4.2) through (4.6) may be written as 

{zO+D(x)} = {z0)(x)} + {AzO)(x)J (4.7a) 

where 


{Az(0(x)}' = ([A(x,z(0)] + [B(x,z(0)]) {Az0)( x )} 

- (y ( u i} (>0) 


(4.7b) 



and 


{y ( u °(x)} = ([A(x,z0-1))] - [A(x,z‘)]) [z0)(x)} 

(4.7c) 

( a . ^ 

( B j} = [A(x,z(i))j (z(0(x)} 

(4.7d) 

(Az0)(x)} = {0} at boundary points 

(4.7e) 


A few iterations of equations (4.7) are required to achieve 
convergence for the nonlinear steady deflections. This procedure is 
essentially equivalent to the quasi-linearization method of solving 
nonlinear boundary value problems [42]. 

4.3 Static Transfer Matrix and Solution 

The starting and the subsequent iteration solutions (equations 
(4.2) - (4.7) involve the solution of a linearized branched boundary 
value problem of the following form. 

{z(x)}’ = [A(x)] {z(x)} + {f} (4.8) 

The transfer matrix for the homogeneous part of this equation is 
static (time independent) and is given by [23] as 

[T(x)]' = [A(x)] [T(x)] (4.9) 

with the initial conditions 


[T(0)] = [I] 


(4.10) 
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The static transfer matrix is used to construct the homogeneous 
solution. Taken together with the particular integral the complete 
nonhomogeneous solution for equation (4.8) can be written as in 
reference [26] 

(z(x)} = [T(x,0)] { z(0) } 

x 

+ [T(x,0)] J [T(s,0)]-1 (f(s))ds (4.11) 

o 

4.4 Formulation 

From equation (4.11), the following relation can be written 
between the state vectors at locations 2 and 3 on the blade (see 
Figure 3.1) 

{ 23 } = [T] {Z 2 } + {c} (4.12) 

This equation may be rewritten in an expanded form relating 
deflections (d) and forces (f). 

T 11 I T 1 2 

| 

T 21 I T 22 

Similarly, the transfer matrix relation for the i th load path can be 
written (see Figure 3.1) 

{4} = [T»] {z}} + {c»} (4.14) 

Expanding the above equation into a partitioned form gives 
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The transfer matrices in the above equations are determined from 
equations (4.9) and (4.10) and the nonhomogeneous solution vector 
{c} is obtained from equation (4.11). 

In general, boundary conditions are defined at the branch roots 
and the blade tip. Equilibrium (equation (3.2)) and compatibility 
(equation (3.4)) relate the state vector between the blade root and 
the branch tips (across the clevis). 

The boundary condition at the blade tip is 

{f3} = {0} (4.16) 

The root ends of the branches are assumed either as clamped in 
bending and extension and either clamped or spring restrained in 
torsion. However, the formulation allows for different root 
conditions. The number of branches is also kept general: 

clamped: {d j } = {0} (4.17) 

spring restrained (torsion): 


m = vi = wi = 0; <> i = M xl /k<J> 
Myi = 0, M z j = 0 


(4.18) 
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— 

where - control system stiffness. 



From equation (4.15) 



{d|} = tT^] {d\) + [T/ 2 ] {f }} + {c\) 

(4.19) 


{4)=[T 2 i 1 ] {d/l + IT^] {fjj + fc!,} 

(4.20) 


From Eq. (4.19) 



{fi) = [T, i 2 ] -*(dj) -[T/j] -Ht l (cj } 

(4.21) 


Apply compatibility across the clevis. Substitute equation 

(3.4) into 

— 

(4.21) 


— 

(fi)=[T, i j-l[B i ] 1 |d 2 ) - [T,' 2 ] _i It /,] {d j ) 



-[T ,2 ] ([B']"’{b) + (cj) 

(4.22) 


Eq. (4.22) for clamped branches can be written by applying 

the root 


boundary condition (4.17) 


— 

{fj} = [r 1 ] {d 2 } + {e i } 

(4.23) 


where 



[r 1 ] = [T/ 2 ] ItB 1 ]- 1 

(4.24) 

— 

{ei)=-[T 1 i ^-l([B i ]- 1 (b} + {cj}) 

(4.25) 


For spring restrained branches, equation (4.22) remains as 


— 

{fj} = [r 1 ] {d 2 } +[s i ] {dj} + {e 1 } 

(4.26) 






Expanding Eq. (4.26) and substituting boundary conditions given by 
Eq. (4.18) yields 
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0 

0 


< 


> = 


V 


v* 


yi 


V 1 

Z1 


V' 

XI 


J 


H l 
... 

H2 

< 


ld 2 1 

1 

_ T21 1 

1 r 22 - 


l 

- - - r 
r d2 


SI 1 
S21 


I si2 
I - - 

I S22 


r o i 


< 


>+< 




Extract the top three equations of Eq. (4.27) and rewrite as 


~ 0 0 
0 0 0 
-000 


R) 


1«1 




= [rill (ld 2 ) + I r 12 ] ( r d 2 } 


Ui 


\ J 


(4.27) 
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+ [S12] 



(4.28) 


Rewriting Eq. (4.28), 


e l 


[p] \ «s 






= [rill {Id 2 } + [ri2] {rd 2 } + {e}} 


\ J 


where 


[p] = 


0 0 k* 

0 0 0 
0 0 0 J 


- [si 2 ] 


From Eqs. (4.29) and (4.18) 


(dj) = [a'] (d 2 ) + {m'j 


where 


[a*] 


[0] 

1 [0] 

. tp]-'[rii] 

1 

1 [pWtai] . 

r. .° . 

.1 


L CpJ-UeJ} J 


(4.29) 


(4. .30) 


(4.31) 


(m*} = 


(4.32) 



Substituting Eq. (4.30) into Eq. (4.26) yields 
{f|} = ([ri] + [si] [ai]) [d 2 ] + [si] [mi] + [ei] 


(4.33) 
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For clamped branches, the boundary condition (4.17) can be applied 


to equation (4.20) 

{f‘)=rry (f,')+ (4) (4.34) 

Substituting Eq. (4.23) into Eq. (4.34) gives 

(4)=[k i ](d 2 ) + (l i ) (4.35) 

where 

[k 1 ] = [Tj 2 ] [r‘] (4.36) 

[1*]= [T 2 ^ [e 1 ] + [c*] (4.37) 


For spring restrained branches, substituting Eq. (4.30) into Eq. (4.20) 
gives 

( 4 ) = [k i ](d 2 ) + (l i ) (4.38) 

[k‘] = [T 2 ‘J [r 1 ] + ([T 2 'j ] + [Tj’jl [s']) + [a'] (4.39) 

[l 1 ] = [Tjjl [e' } + [cjj + dTj*,] + [T' 22 | [s 1 ] ) (m 1 ) (4.40) 

Substituting Eq. (4.35) or Eq. (4.38) into the equilibrium equation 
(3.2) across the clevis yields 

[h] = [ki] [d2] + [k2] (4.41) 


where 
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n 

[kl] = X [A 1 ] [k'l ( 4 . 42 ) 

i=l 

n 

M= ItAi] ( ii) (4 43) 

i=l 

From Eq. (4.13) 

{ f 3 } = U21] {d2} + [T22] {f2} + { C2 } ( 4 . 44 ) 

Substituting Eq. (4.41) into Eq. (4.44) yields 

{ f 3 > = ([T21] + [T 22 ] [ki]) {d 2 } + [T 22 ] {k 2 } + (c 2 ) ( 4 . 45 ) 

From applying the boundary condition for the blade (Eqs. (4.16) - 
(4.45)) 

(d 2 ) = -([T21] + [T22] [kj ])' 1 ([T22] {k2} + {c2 ) ) ( 4 . 46 ) 

Now the state vector [ 22 } is known from Eqs. (4.46) and (4.41). The 
state vector {zj} for the spring restrained branches can be obtained 
from Eqs. (4.30) and (4.33). For clamped branches {dj} = (0) by 
virtue of the boundary condition and {fj} can be computed from Eq. 
(4.23). Once the state vectors {zj} and {Z 2 } are known, the solutions 
at any station can be determined from Eq. (4.11) as shown below. 


i th branch: 

{z‘(x)} = [T*(x)] {zj } + {c'} 


(4.47) 
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blade: 

{z(x)}=[T(x)} i^ + ic) 


(4.48) 


For blades with single root branches, the boundary conditions are 
assumed as shown below. 


x = 0: u = w = v = e = ^= 0 
x = R: {f 3 } = {0} 



(4.49) 

(4.50) 


The equilibrium and compatibility matrices become identity matrices 
and this yields 


(d 2 ) = {dl} = [T /,] {dj ) + [T > 2 ] {f, } + {c,} (4.51) 

(f 2 )= (f|) = lT 2 ' 1 ] Id,) + [T 2 2I (f,) + {03} (4.52) 

Substitute Eqs. (4.50) to (4.52) into Eq. (4.45) yields 

[A]{di)+[B](f,) + (E) = (0) (4.53) 

where 


[A] = [T 2i ] [T 1,1 + [t 22 ] [T 22 ] 

[B] = [T 21 ] [T / 2 ] + [T 22 ] [T 2 * 2 ] 

[C] = [T 21 ][cj] + [T 22 ] [c 2 ] 


Substituting Eq. (4.49) into Eq. (4.53) yields the following equation 
for 4>i. 
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<1>1 = -La 0 0 0 0 Oj[B]-l{E) (4.54) 

where 

a = 1 /(Di ,6 + k<f>) 

and 

[D] = [B]-l [A] 

Once 01 is known, {di } is known and {fi } can be computed from Eq. 
(4.53) as shown below. 

{fl} = -[B]-l[A]{di)-[B]-l[E] (4.55) 

Once the state vector { zi } is known, the solution at any station can be 
computed from the following equations. 

branch: 

{ z(x) } = [T*(x)] (zi) + (c) (4.56) 


blade: 

(z(x)} = [T(x)] { z 2 } + {c} (4.57) 

4.5 Solution Procedure 

As described in section 4.2, the linear solution is taken as the 
starting solution for the iterations. The linear solution is obtained by 
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solving Eq. (4.1) with (z(x)} = {0} in the coefficient matrix except for 
the tension V x . The tension in the blade is obtained as shown below. 

R 

V x (x) = Jn 2 mxdx (4.58) 

x 


The tensions in the branches are calculated by assuming that the 
branches are coincident with the blade at the clevis, i.e., hj = hj, = 0. 

The tensions corresponding to this case are calculated as follows. The 
transfer matrix for axial motion of a beam for the static case is given 
by 


1 


[T(x)J = 


a 

1 


(4.59) 


where 



By definition of the transfer matrix 



(4.60) 


From the above equation 



(4.61) 
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The following two cases are considered 

Case 1: 2 branches 

Expansion of Eq. (4.61) yields 


Uo = 


u? = 


a I V x 2 

a 2 V * 2 2 


For coincident nodes (u^ = u^), Eq. (4.62) becomes 


a 1 V x 2 - a I V x 2 = 0 


For equilibrium 


V 1 + V 2 = v 

x 2 x2 V X2 


Solving Eqs. (4.63) and (4.64) gives 


V 1 = — — V 

x 2 ai + a 2 x 2 

V 2 = — — V 

x 2 ai + a2 x 2 ) 


Case 2: 3 branches 

Expanding Eq. (4.61) for coincident nodes yields 


a l V J2 =a 2 V x 2 2= a 3 V x 3 


x 2 


(4.62) 


(4.63) 


(4.64) 


(4.65) 


(4.66) 


For equilibrium 
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V 1 + v 2 + V 3 = v 

x 2 x 2 v x 3 v x2 
By solving Eqs. (4.66) and (4.67) yields 


^X2 a 2 a 3^x2^^ 

V x2= a 3 a , V x 2 /D 

V x2 =a l a 2 V x 2 /D, 


where D = aia 2 + a 2 a 3 + a 3 aj 

Now the above result can be generalized to the ' 
as shown below. 


V 


i 

x 2 


n,i 

X a j 

hi 

n n,k 

X <X aj) 

k= 1 j=l 


where 



a n/ a i 


and V X2 is obtained from Eq. (4.58). 

Now the tension in the ith branch corresponding to the 
branch case is given by 

1 

Vic0 = J Q 2 mxdx + V^ 2 

X 


(4.67) 


(4.68) 


branch case 


(4.69) 


coincident 


(4.70) 



Once V x in the branches and the blade is known, the linear steady 
state solution can be obtained following the procedure outlined in 
Section 4.4. Note that for the single branch case, the V x is simply 
given by equation (4.58) in both the branch and blade. 
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5. LINEAR PERTURBATION EQUATIONS 

In the absence of experimental flight testing, the aeroelastic 
characteristics of any flight vehicle in trimmed flight may be 
investigated using perturbation equations. The procedure for the 
development of perturbation equations is universal for all flight 
dynamics problems. Here the complete perturbation equations for 
an investigation of the flight dynamic instabilities of an isolated rotor 
blade are presented. Components of such an aeroelastic analysis are 
i) free vibration characteristics about the nonlinearly deformed 
steady state and iii) complex stability eigenvalues. In the present 
work, the equations below are utilized to investigate each of these 
three elements. 

The procedure for the development of the linear perturbation 
equations is 

1. Substitute (z(x,t)} = (zo(x)} + (z p (x,t)} into the governing 
differential equations (2.7), (2.8), (2.9) and (2.33-2.41) 
and 2.80-2.82. 

2. Subtract the nonlinear steady state equations (4.1) from 
the results of step 1. This step will eliminate all steady 
state and nonhomogeneous terms. 

3. Drop the nonlinear terms in the perturbed variables 
{z p (x,t)}. 

The resulting linearized homogeneous perturbation equations are 
listed below (with the bars omitted for simplicity). 
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u' = -e 0 e - + V X /EA 

w' = e 
v' =% 

e' = -aiM Xo e + a 3 M X( ^ - (2aiMy 0 + a 4 M Zo )<(> 
+ (-aie 0 + ai^ 0 )M x + (ai - a4<t>o)M z 
- (a3 + 2ai<t> 0 )My 


^ — -a 2 M XQ e + aiM Xo ^ — (2aiM Zo - a 4 My o )<]) 
+ (-a 2 e 0 + ai^ 0 )M x 
+ (a 2 - 2ai<t> 0 )M z - (aj - a 4 <t> 0 )My 


<j)' = (1/GJ) (M Zo e + M y ^ + M x + e 0 M z + ^ Q M y ) 

M x = -2meQsin0u + me(cos0-<t> o sin0(w 
- me(sin+<j> o cos0(v - Q 2 v) 

+ 2Qm(k^ 2 sin0 + k^cos^e + V yo £ 

+ 2Qm(k 2 2 ~k 2 ^cosOsinO^ - V Zq $ + mk 2 '<j> 

+ £2 2 m(k^ 2 - k m 2 j) cos 2 0 - <}>osin 2 0) + me£l 2 v o cos0}(|) 

+ eoVy - S 0 V Z - Me 


(5.1) 

(5.2) 

(5.3) 

(5.4) 

(5.5) 

(5.6) 

(5.7) 


M z = 2meQcos0v + V Xq 2; - meQ 2 xsin0<|> - V y + £ 0 V X 


(5.8) 



66 


My = 2meQsin0v - V Xo e - me£2 2 xcos0<j> + V z - e 0 V x (5.9) 

Vy= 2Qmii + m(v-£22 v ) - 2f2mesin0e 

- 2Qmecos04 + mesinOCQ^.^) - L v - 2Qm(3pcw (5.10) 

Vz = mw + mecos0<j> - L w + 2mQ|3p C v (5.11) 

Vx = m(u-f2 2 u) - 2i2mv (5.12) 

where the perturbed aerodynamic forcing functions are 

L v = k 4 [-QxVj o <|> - {2Qxk2+ (9+<J>o) v i 0 ) v 
+ {2v io - £2x(0+<J) o ) io } w 

- kiv io (<{)+f2e) - aki(e 0 +Ppc) io }w] (5.13) 

Lw = k 4 [Q 2 x 2 4 > - fl 2 x(v 0 e+e 0 v) - Q^x(3 pc )v 

+ kiii 2 xe + {2Qx(0+<)>o) - v io + Qki(e<)+Ppc)}v 

- Qx w + kiQx(f>) (5.14) 

M e = k3L w - k 4 C 2 [Qx }<^+ Q 2 xe + Q(eo+P pc ))}v (5.15) 


and 


ki = c(l-e)/2 
k2 = cd/a 


k 3 = ce/2 
k 4 = pac/2 
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6. FREE VIBRATION CHARACTERISTICS 

6.1 Equations for Free Vibration 

The free vibration equations are a special case of the complete 
linear perturbation equations developed in Chapter 5 (equations 
(5.1) - (5.12)). The free vibration equations are developed from the 
linear perturbation equations by the following procedure. 

1 . Drop damping type terms, i.e. terms containing first 
derivatives with respect to time (u,v, w, e, <J>) 

2. Eliminate time dependency in the equations by assuming 
simple harmonic motion with frequency « for u, w, v and 
<{>. 

(z(x,t)} = Ze Ib)t 

3. Drop all the aerodynamic loadings (forcing functions). 

The resulting equations are summarized below: 


u’ = -e 0 e - ^ + Vx/EA 

(6.1) 

w’ = 8 

(6.2) 

v' = 5 

(6.3) 

e' = aiM XQ £ + a 3 M XQ ^ - ( 2 aiM yQ + a 4 M ZQ )<j> 


+ (-aie 0 + a 3 ^ 0 )M x + (ai - a4<}>o)M z 


- (33 + 2ai<J)o)My 

(6.4) 
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5’ = - a 2 M Xo e + aiM Xo £ + (a 4 M yQ - 2aiM ZQ )<|> 

+ (a l^o - a2E 0 )M x + (a2 - 2ai<|>o)M z 
+ (~ai + a 4 <|) 0 )My (6.5) 

= (l/GJ)(M Zo e + My o ^ + M x + e 0 M z + eoM y ) (6.6) 

M x = -co 2 me£2(cos0 - <t> 0 sin0)w + me(sin0 

+ <)>o cos0)(co 2 + £2 2 )v + V v e - V- £ 

jo 

+ {-0) 2 k^ + ^ 2 m(k^ 2 -k^ j )(cos20 - <|)osin20) 

+ meQ 2 v o cos0}<|> + e 0 v y - £ 0 v z (6.7) 

M; = V Xq £ - meQ 2 xsin0(|) - V y + £ 0 V X (6.8) 

= -V XQ e - meQ 2 xcos0<j) + V z - e 0 V x (6.9) 

V y = -(co 2 + Q 2 )mv + (g) 2 + Q 2 )mesin00 (6.10) 

V' z = -G) 2 mw - meG) 2 cos0(j) (6.11) 

V x = -(G) 2 + Q 2 )mu (6.12) 


Equations (6.1) to (6.12) govern the free-vibration state about the 
nonlinear steady state position. Equations of motion governing the 



free-vibration state about the initial geometry can be obtained by 
substituting the following relations in the above equations. 


U 0 = Vo = w 0 = e G = = <}>o = M = M v = M 7 =0 

A o / o 


(6.13) 


R 


V = [ Q 2 mxdx 

* 0 J 

x 

Vy o = Q 2 mxecos<l> 
V z = Omxesin|3 




(6.14) 


6.2 Dynamic Transfer Matrix 

The first-order differential equations of motion are linear and 
homogeneous and therefore can be arranged into a matrix differ- 
ential equation of the following form: 

{z(x)}’ = [A(x)] [T(x>] (6.15) 

The transfer matrix for the above system is given by solving the 
following equations. 

[T(x)]' = [A(x)] [T(x>] (6.16) 

[T(0)] = [I] (6.17) 


Once the transfer matrix is known the free-vibration characteristics 
of branched blades can be determined in the following section. 
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6.3 Frequency Determinant 

To compute the natural frequencies of the branched blade, a 
frequency determinant is computed and scanned for different input 
frequencies until the value of the determinant is zero. The input 
frequency that returns a zero value for the frequency determinant is 
a natural frequency of the system. Since the branches and blade are 
modeled separately, the frequency determinant is constructed by 
relating the state vectors at the branch roots to the state vector at 
the blade tip. 

By definition of the transfer matrix, the following relation 
between state vectors at locations 2 and 3 can be written (see Figure 
3.1) 


{Z 3 ) = m {z 2 } (6.18) 

Rewriting this equation into the following partitioned form: 


d3 

< - _ _ 


. f3 . 



Til I T12 

| 

T21 I T22. 



' d 2 


. f2 . 


(6.19) 


Extracting the following equation for forces from Eq. (6.18) yields 

{f 3 } = [T 2 i] {d 2 } + [T 22 ] (f2> (6.20) 

Similarly, the transfer matrix relation for the ith branch can be 
written as (see Figure 3.1) 
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(4) = [H] (zj) 

Rewriting the above equation into partitioned form: 



( 6 . 21 ) 


( 6 . 22 ) 


Expanding Eq. (6.22) gives 

{d*} = [T 1 i 1 ] {d*} +[T, i 2 ] {f>) (6.23) 

{ 4 ) = [T 2 i i) (dj) + [Tj 2 ] (f{) (6.24) 


The boundary conditions for the branch roots are given by 

Clamped: {di } = 0 (6.25) 

Spring Restrained: ui = vj = wi = M yi = M Z] = 0 (6.26) 

For homogeneous problems, the collective pitch can be added to the 
pretwist of the blade and this will translate as elastic twist to the 
branches. For free-vibration problems, the steady state elastic twists 
of the branches have to be determined to solve the homogeneous 
problem. Once this is done, the compatibility relation given by Eq. 
(3.4) can be written as 


{d 2 } = [Bi] {d‘ } 


(6.27) 



If the free-vibration is solved about the undeformed state, the elastic 
twist of the branches due to collective pitch should be included as 
pretwist of the branches. The equilibrium equation across the clevis 
will be same as before as given by Eq. (3.4). 

n 

<f 2 ) = I [A*] (4) 
i=l 

Two cases are considered here (clamped and spring restrained 
branches). 


1 . Clamped branches 

By virtue of the boundary condition equation (6.25), equations 
(6.23) and (6.24) become 

(4)=IT 1 i 2 ] (fj) (6.28) 

f f 2)=[ T 22] f f j) (6.29) 

From Eq. (6.28) 

{fj) = [T jj' 1 (dj) (6.30) 


Substituting the compatibility Eq. (6.27) into Eq. (6.30) gives 
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Substituting Eq. (6.31) into Eq. (6.29) yields an expression for forces 
in the clamped branches at the clevis 

{ 4 } = [T 2 ‘ 2 ] [Tj'j-l |B]-‘ (d 2 ) (6.32) 

or (f‘) = [k'] (d 2 ) 

where [k 1 ] = [T 2 ‘ 2 ] [r 1 ] (6.33) 

[r‘] = [T |j l [B*]- 1 (6.34) 

2. Spring restrained branches 
From Eq. (6.23) 

(fi)=( T , i 2 ) -Mdj) -[Tj' 2 ]-' [ T ,>,] (dj) (6.35) 

Substituting the compatibility Eq. (6.27) into Eq. (6.35) gives 

{f»] =[T*^ -1 [B 1 ]- 1 {d 2 } - [T,' 2 ]-> [ T\ ,] (dj) (6.36) 

Rewriting Eq. (6.36) as 

Ifi)=[r i J{d 2 ) + ts*] (dj) 


(6.37) 
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where 

[r 1 ] = [Tjj-1 

[s 1 ] =-[T{j-l [ T },] 


Expanding Eq. (6.39) 




K =v»i 
0 
0 


V 1 
yi 


V 1 

zi 


> - 


V! 


XI 


ri 1 I ri2 


‘d2 


L T21 | T22 j 


r d2 J 


si 1 | Si 2 
- - I - - 

L S21 | S22 j 


0 

0 

0 


*i 


(6.38) 

(6.39) 


Extracting and rewriting the top three equations of the above 
equation yields 
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■ 0 

0 

k < j ) 


fVl 

E i 

0 

0 

0 

< 

Si > 

_ 0 

0 

0 _ 


<*\j 


*d 2 J 


r d 2 


r i ^ 
e l 


+ [S12] 


UJ 




From Eq. (6.40) 


r 






> _ 


= [p]-> [ru] (l d2 ) + tpJ 1 [T12] (r d2 ) 


-1 



where [p] is defined in Eq. (4.29) 

From Eqs. (6.41) and the boundary condition equation (6.26) 
spring restrained branches 

(dj) = [a‘] (d 2 ) 


where 


[ai] = I 


0 


L [p] _1 [ri i ] | [p]" 1 [ri 2 ] J 


(6.40) 

(6.41) 
for 

(6.42) 


Substituting Eq. (6.42) into Eq. (6.41) yields 
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(fi)=[r i ]{d 2 ) + [s i ][a i nd 2 } (6.43) 

Substituting Eqs. (6.42) and (6.43) into Eq. (6.24) yields 

(4) = [k 1 ] {d 2 } (6.44) 

where 

Ck 1 ] = [T^ [r 1 ] + ([Ty [s 1 ] + [T^j]) [ a i] (6.45) 

Equilibrium is enforced by substituting Eq. (6.44) into Eq. (3.4), 
giving 

{f2} = [D]{d 2 } (6.46) 

where 

n 

[D]= X [Ai] [ki] (6.47) 

i=l 

Substituting Eq. (6.46) into Eq. (6.20) yields 

(f3) = ([T 2 i] + [T 22 ] [D]) {d 2 } (6.48) 

The boundary condition for the blade tip is given by 

{f3} = {0} (6.49) 


Application of this boundary condition to equation (6.48) yields 
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([T 2 1 ] + [T 22 ] [D]) {d 2 } = {0} (6.50) 

For nontrivial solutions of {d 2 }, the determinant of the coefficient 
matrix should be zero and this condition yields the following 
frequency equation to determine the natural frequencies of the 
multiple branch blades. 

det([T 2 i] + [T 22 ][D]) = 0 (6.51) 

For single branch blades, the boundary conditions are assumed as 


shown below. 

x = 0: ui = wi = vi = ei = £i = 0 (6.52) 

<t> l = M xl /k ( | ) 

x = R: {f 3 } = {0> (6.53) 

The equilibrium and compatibility matrices become identity matrices 
and this yields 

{d 2 ! = (d 2 ) = [Tj,] (d,) + [Tjj] (fj) (6.54) 

(f 2 ) = (f 2 ) = [T!,] (d, ) + [ T> 2 ] (f,) (6.55) 

Substituting Eqs. (6.54) and (6.55) into Eq. (6.20) yields 

{f 3 } = [A] {d^ + tB] (fj) (6.56) 


(6.52) 
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where 

[A] = [T 21 ] [T | j] + [T 22 ] [T-, 1 , ] 

[B] = [T 21 ][T, 1 2 ] + [T 22 ][T 2 1 2 ] 

Application of the boundary condition equations (6.52) and 
equation (6.56) gives 


r o > 


r M x = M' 

0 


M z 

0 


Mv 

< 

0 

” + [B] < 

y 


Vy 

0 


v z 

^<{» 

1 

L y x J 


From Eq. (6.59) 


[D] (fi) = 


- < 


r 0 ' 


0 

0 

0 

0 


♦l 


(6.57) 

(6.58) 

(6.53) to 

(6.59) 


(6.60) 


where 



Extracting the last row of the above equation yields 


D 6,l k <|)^l + L D 6,2 D 6,3 D 6,4 D 6,5 D 6,6 ^ 


/m z \ 

My 

Vy 
V z 
VVx^l 


= <t>i 


Rearrange Eq. (6.61) as 


<h i — L a 4 a 5 a g J { f | } 

where 

D,. 

a = 0 and a. = — 0 : , j = 2 to 6 

1 J D 6.lV _1 

Substituting Eq. (6.42) into Eq. (6.60) yields 
([D] - [A]) {fi} = {0} 


(6.61) 


(6.62) 


(6.63) 


where 
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For nontrivial solution of {fj}, the determinant of the coefficient 

matrix of Eq. (6.63) should be zero and this condition yields the 
following frequency equation to determine the natural frequenices of 
the single branch blade. 

det([D] - [A]) = 0 (6.64) 


6.4 Mode Shapes 

A mode shape may be defined as the shape corresponding to a 
specific frequency in which the elastic and inertial forces are in 
equilibrium. The formulation for calculating the fully coupled flap, 
lag and pitch mode shapes follows. 



Til I T 12 


T21 I T22 



where 


(6.65) 


[T] = [T]-l 
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From the above equation, the displacement vector at the clevis is as 
shown below by virtue of the boundary condition at the blade tip 

(f 3 ) = {0}, 

{d 2 } =[T n J {d 3 } (6.66) 

Substituting Eq. (6.66) into Eq. (6.48) yields 

{f 3 } = [C]{d 3 ) (6.67) 

where 

[C] = ([T U ] + [T 22 ] [D]) [T n ] 

Assuming w 3 = 1 arbitrarily and rewriting rows 2 to 6 of Eq. (6.67) 
gives 


C 2,l C 2,3 02,4 C 2,5 
C3,l C 3>3 C 3 ,4 C 3> 5 
Q,1 C4,3 C4.4 C4,5 
C5.1 C5, 3 C5,4 Cs,5 
C6.1 C6.3 C6,4 C6.5 


c 2 , 6 " 


u 3 


r -c 2 , 2 " 

C 3 ,6 


v 3 


-C3, 2 

C4.6 

< 

e 3 


-c 4>2 >• 

C5.6 


^3 


-C5,2 

C6.6 - 




L--C6, 2 J 


( 6 . 68 ) 


By solving the above equation, u 3 ,v 3 ,e 3 ,^ 3 , and <|> 3 are known and 
together with w 3 = 1 the entire blade tip deflection {d 3 } is known. 
Once {d 3 } is known, the state vector at the clevis can be determined 

ferom Eq. (6.65) as shown below. 



Once the state vector at the clevis is known, the deflection vectors in 
the blade and the branches can be obtained as follows. 


Blade: By definition of the transfer matrix, the state vector at any 

location x is given by 


' d x 
- 

> = 

"Tn 

T 2 f 

— 


' d2 ■ 

. fx . 


_T 2 1 

1 T22_ 

X 

. f2 , 


(6.70) 


From the above equation 


( d x)=[Ti,] x {d 2 ) + [T,2] x (f 2 ) 


(6.71) 


Branches: By definition of the transfer matrix, the state vector at 

any location x in the branch is given by 


J “M 

‘ T A 

1 T 12 


d i 1 

— f = 

- - - 

1 - - - 

< 

- - - r 

I f i J 

- T ii 

i n 2 ] 


^ \ j 


For clamped root branches the following equation can be extracted 

from Eq. (6.72) by virtue of the boundary condition deflections 
{dj } = {0} at the root. 



83 


{d*J = f T i i 2 ] ( f i> (6-73) 

For spring restrained root branches the equation remains as 

{dl) = [T l \) {dj} + [T/ 2 ] {fj } (6.74) 

Eq. (6.31) provides {fj} for clamped branches, and Eqs. (6.42) and 
(6.43) provide {dj),{fj} for spring restrained branches. 

For single branch blades Eq. (6.63) can be solved for { f i } 
assuming V z = 1 arbitrarily, and then <)>i can be comnputed from Eq. 
(6.61). 

branch: (d x ) = [Tj, ] (d,) + [T 1 ^] (f,) (6.75) 


blade: 


(d x )=[T„] {d 2 )+ [T 12 ] (f 2 ) 


(6.76) 
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7. AEROELASTIC STABILITY CHARACTERISTICS 

7.1 Equations for Aeroelastic Stability 

The perturbation equations (5.1) - (5.12) are specialized with a 
transformation to the frequency domain to determine the equations 
governing the aeroelastic stability characteristics of branched blades. 
The transformation is accomplished by substituting 

{zP(x,t)} = ZP(x)} e*t (7.1) 

where X is complex and equal to o + ico . (Note that for free vibration X 
was simply equal to ico.) The resulting equations are (superscripts p 
omitted) 

u' = -e 0 e - ^ + V X /EA 
w' = e 

v’M 

e ’ = -aiM Xo e + aiM x ^ - (2aiM yQ + mM Z o )$ 

+ (~aie 0 + a3^ 0 )M x + (ai - a4<])o)M z 

- (a3 + 2ai<()o)My (7.5) 

% = * a 2 M x 0 e + a l M x 0 ^ - < 2a l M z 0 “ 34M yo )(t> 

+ (— a2£o + ai^ 0 )M x + (a2 - 2a]<J>o)M z 

- (ai - a4<))o)My (7.6) 


(7.2) 

(7.3) 

(7.4) 


o-a . 
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<t>' = (l/GJ)(M Zo e + M y ^ + M x + e 0 M z + ^M y (7.7) 

M x = -2me£2sin0^u + me(cos0-<])o sin0)X, 2 w 

- me(sin0+<J)o cos0)^ 2 -£2 2 )v 

+ 2Qm(k^ 2 sin 2 0 + iCOS 2 0 )^ e + y^ e 

+ 2Qm(k^ 2 -k^ } )cos0sin0X^ - V Z J, 

+ {mk^ 2 X 2 + Q2 m (k2 ^k^Xcc^O - <|>osin20) 

+ meQ 2 v o cos0}<j> + eoVy - ^ 0 V Z - Me (7.8) 

M^ = 2QmeQcos0^v + V Xo £ - meQ 2 xsin0<J> - V y + £ 0 V X (7.9) 

M^ = 2Qme£2sin0Xv - V Xq £ - me£2 2 xcos0<|> + V z - e 0 V x (7.10) 

Vy = 2QXu + m(A, 2 -Q 2 )v - 2Qme^.sin0e 

- 2QmeXcos0^ + mesin0(Q 2 -A, 2 )<|> 

-L v -2mQppc^w (7.11) 

V^, = mX 2 w + meX 2 cos0(|) - L w + 2mQp pc v (7.12) 

V x = m(k 2 -Q 2 )u - 2Qm^.v (7.13) 


where 
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L v = k4[-Qxvj o <J) - {2Qxk2 + (e+<t>o)Vj }^v 
+ {2 vj o -Qx(0+<|> o )Xw 


- klv io (X<t>+f2e) - Qki(e 0 + P pc )kw 

(7.14) 

= k 4 [£ 2 2 x 2 <J) - Q 2 x(v 0 e+e 0 v) - £2 2 xp pc v 


+ kiQ 2 xe + {2Qx(6+4> 0 ) - v iQ + Qki(e 0 +P pc )}?tv 


- QxA,w + klQxtaj)] 

(7.15) 

= ^ 3 L w - k4C 2 [f2xX.<t> + Q 2 xe + QA,(e 0 +Ppc) v l 

(7.16) 


and 

kj = c(l-e)/2 k 3 = ce/2 

^2 = Cd/a k 4 = pac/2 

7.2 Stability Eigenvalues 

Equations (7.2)-(7.13) are linear homogeneous equations like 
the free vibration equations (6. 1 )-(6. 1 2), and calculation of the 
complex transfer matrices proceeeds similarly for the branches and 
blade. Equations (7.2)-(7.13) can be arranged into a matrix 
differential equation 

{z(x)}’ = [A(x)]{z(x)J (7.14) 

where [A(x)] is now complex. The transfer matrix for the above 
system is given by solving the following equations 
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[T(x)]' = [A(x)] [T(x)] (7.15) 

[T(0)] = [I] (7.16) 

Once the transfer matrices are known for the blade and branches a 
stability determinant is formulated. The stability determinant is of 
the same form as the frequency determinant (equation 6.47) and is 
given by 

det([T 2 i] + [T 22 ] [D]) = 0 (7.17) 

where [D] is defined with equation (6.43). 

The eigenvalues of the stability determinant are complex 
whereas the eigenvalues of the frequency determinant are real. 

Thus X in equations (7.2) - (7.16) takes complex values. The 
eigenvalues of the stability determinant are calculated by Muller's 
method [42] for determining complex roots, whereas the eigenvalues 
of the frequency determinant are calculated by a frequency scanning 
technique. 

Stability of the branched blade is inferred by examining the 
sign of o (the real part of X). If a is positive then it can be seen from 
(7.1) that the perturbation state vector increases exponentially with 
time and thus the system is unstable. If a is negative, then the 
perturbation state vector decreases with time and eventually damps 
out. In this case the system is stable. 
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8. NUMERICAL RESULTS AND DISCUSSION 

8.1 Computer Program 

A comprehensive computer program is developed to imple- 
ment the formulation presented in the preceding chapters. The 
program determines the following for a branched blade with up to 
three branches: 

1. Natural frequencies and mode shapes about the 
undeformed position 

2. Steady state nonlinear deflections of the blade and 
branches 

3. Natural frequencies and mode shapes about the steady 
state deformed position 

4. Complex stability eigenvalues 

Additionally, the program features fully coupled nonlinear 
flapwise bending, chordwise bending, torsion and axial extensions. It 
can also handle noncoincident mass, elastic and aerodynamic center 
axes with nonuniform property distributions in both the blade and 
branches. The code is lengthy (approximately 4000 lines) and thus is 
included in a separate volume. 

The program is based on a continuous system model, with 
transfer matrices calculated by a fourth order Runge-Kutta 
integration scheme. It should be noted that if a discrete model is 
used to compute the transfer matrices the formulation is still valid. 
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The transfer matrices based on the discrete model can be used in 
place of those computed based on the continuous system model. 

8.2 Validation of Nonlinear Formulation - Single Branch 

Blade 

To demonstrate the extension of the transfer matrix method to 
nonlinear problems, the conventional single-branch rotor blade 
model considered in reference [40] is analyzed first. The data of the 
model is shown in Table 8.1. 

The natural frequencies (both rotating and nonrotating) are 
given in Table 8.2. Nonrotating frequencies about the initial 
undeformed state are computed using the present method and closed 
form analytical expressions (see Appendix B). The rotating blade 
frequencies are shown about both the initial undeformed position 
and the nonlinear steady deformed position. The flap frequencies 
are strongly effected by rotation due to added stiffness coming from 
centrifugal effects. Chord and torsion frequencies are much less 
sensitive to rotation. The frequencies about the deformed state are 
seen to be close to those about the initial undeformed state. 

The steady state nonlinear tip deflections are given in Table 
8.3. These trim tip deflections are compared graphically with the 
results of reference [39] in Figure 8.1, and the agreement is quite 
good. The Newton-Raphson iteration scheme developed for the 
nonlinear distributed system equations is employed to obtain the 
convergence. The efficiency of the scheme is indicated in Table 8.4 
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Table 8.1 

Data for conventional blade 


C/R = 40 in 
0 = 0 . 0 ° 

P pc = -2.8648o to 5.72960 
cow = 1.15 
CDv = 1.50 
city = 5.0 

y = (3paCR/m) = 5 


km2 — 0 
k m2 /2 = 0.025 
(k A /k m ) 2 = 1-5 

a = 2 jt per rad 
b = 4 

Cd 0 /a = 0.01/27C 
6 rc = 0.0° to 28.6479° 


e/R = 0.0 
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Table 8.2 

Natural frequencies (rad/sec), single branch blade 


Mode 

Initial state 


Deformed state 


Q = 

0 Q 2 

= 9.3456 

Q 2 = 9.4356 


Present 

Method 

Analytical 

P rc = 0.3 rad 

1 

3.3813 F 

3.3810 

9.1576 

8.5005 

2 

11.4324 C 

11.4314 

11.9370 

12.0323 

3 

21.1905 F 

21.1898 

29.3043 

28.9895 

4 

37.6596 T 

37.6601 

39.7894 

40.3584 

5 

59.3348 F 

59.3381 

68.1363 

67.5846 

6 

71.6468 C 

71.6443 

74.0281 

73.8076 

7 

112.9801 T 

112.9803 

116.4039 

1 17.291 1 

8 

116.2750 F 

116.2725 

125.6590 

124.7716 

9 

188.3004 T 

188.3004 

193.6040 

194.0339 

10 

192.1118 F 

192.2055 

201.9584 

200.5757 


F = Predominantly flapwise 
C = Predominantly chordwise 
T = Predominantly torsion 
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Precone 
Ppc (rad) 

-0.05 (-2.86°) 
0.00 
0.05 
0.05 


Table 8.3 

Steady state deflections, conventional blade 

Root Collective w 0tip v 0tip 4>o tip 

0 rc (rad) 


0.0 

0.1 (5.73°) 

0.2 

0.3 (17.19°) 
0.4 

0.5 (28.65°) 

0.0 

0.1 

0.2 

0.3 

0.4 

0.5 

0.0 

0.1 

0.2 

0.3 

0.4 

0.5 

0.0 

0.1 

0.2 

0.3 

0.4 

0.5 


0.4079x10"! 

0.6197x10'! 

0.9527x10-! 

0.1329 

0.1721 

0.2096 

0.0000 

0.2150x10'! 

0.5554x10'! 

0.9466x10'! 

0.1366 

0.1790 

-0.4079x10-1 

-0.1935x10-1 

0.1503x10-1 

0.5516x10*1 

0.9922x10-1 

0.1454 

-0.8157x10-1 

-0.6058x10-1 

-0.2628x10-1 

0.1440x10-1 

0.6005x10-1 

0.1094 


-0.5266xl0- 3 

-0.6360xl0- 2 

-0.1944x10-1 

-0.4109x10-! 

-0.7189x10-! 

- 0.1110 

-0.5222x10-3 

-0.2942x10-2 

-0.1272x10-! 

-0.3145x10-! 

-0.6024x10-! 

0.9938x10'! 

-0.5266x10-3 

0.5137x10-3 

-0.5813x10-2 

-0.2128x10-! 

-0.4746x10-! 

-0.8553x10-1 

-0.5398x10-3 

0.3977x10-2 

0.1234x10-2 

-0.1069x10-1 

-0.3375x10-1 

-0.6994x10*1 


0.2027x10-3 

-0.4760x10-2 

-0.9353x10-2 

-0.1391x10-1 

-0.2001x10-1 

-0.3180x10-! 

0.0000 

-0.4693x10"2 

-0.8666x10-2 

-0.1 174x10-1 
-0.1469x10-1 
-0.2022x10-1 

0.2027x10-3 

-0.5085x10-2 

-0.8956x10-2 

-0.1120x10-1 

-0.1202x10-1 

-0.1313x10-1 

-0.4070x10-3 

-0.5938x10-2 

-0.1022x10-1 

-0.1224x10-1 

-0.1159x10-1 

-0.9264x10-2 


®olttp 




Table 8.4 

Convergence of nonlinear steady state trim (tip deflections) 
P pc = 0.05 (2.86o), 6 rc = 0.3 rad (17. 19°) 
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State 

variable 

Starting 

solution 

zO)(x) 

Starting 

solution 

z(2)(x) 

I 

iteration 

II 

iteration 

III 

iteration 

— 

u 

0. 3946xl0- 6 

-0. 1913xl0- 2 

-0.2023xl0- 2 

-0.2007x10-2 

-0.2007x10-2 


w 

0.5 140x10"! 

0.5563x10-! 

0.5516x10-1 

0.5516x10-1 

0.5516x10-1 


V 

-0. 2037x10-! 

-0.2145x10-1 

-0.2128x10-1 

-0.2128x10-1 

-0.2128x10*1 

— 

e 

0.7444x10-! 

0.7972x10-! 

0.7912x10-1 

0.7912x10-1 

0.7912x10-1 



-0.2868x10-! 

-0.2990x10-! 

-0.2967x10-1 

-0.2967x10-1 

-0.2967x10-1 


4> 

-0.1406x10-! 

-0.1034x10*1 

-0.1121x10-1 

-0.1120x10-1 

-0.1120x10-1 


95 


which shows the nonzero elements (deflections) of {z} at the blade 
tip. The first and second columns of this table are given by 
equations (4.1a) and (4.1b) respectively. The third, fourth and fifth 
columns are the first, second and third Newton-Raphson iterative 
solutions respectively computed using equations (4.7). All the blade 
tip motions are converged to four significant figures. 

The complex stability eigenvalues are determined from 
perturbations about the nonlinearly deformed steady state position. 
The stability eigenvalues obtained using the present method are 
compared with the results obtained in references [40] and [44] in 
Table 8.5. Since the real parts of the eigenvalues are negative for all 
the modes (flap, lag and torsion) the perturbation motion damps out 
and the system is stable. The flap and torsional motions have a 
much larger stability margin than the lag motion since the real parts 
of those eigenvalues are much larger negative numbers than the real 
part of the lag stability eigenvalue. It is also interesting to note that 
the imaginary parts of the stability eigenvalues in Table 8.5 are quite 
close in value to the natural frequencies computed about the non- 
linearly deformed steady deformed position (Table 8.2). Since the 
imaginary part of the stability eigenvalue represents the damped 
natural frequency and it is quite close to the undamped natural 
frequency it is clear that there is very little damping in this system. 
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Table 8.5 

Comparison of stability eigenvalues, conventional blade 


9 rc = 0.3 rad (17.19°), p pc = 0.0 rad 

Present TM Ref. 40 Ref. 44 

Formulation 
(rad/sec) 


Flap -2.2509 + 7.3895i -2.2449 + 7.5554i -2.2194 + 7.5957i 

Lead-lag -0.4851 + 12.5397i -0.5029 + 12.5915i -0.5308 + 12.4642i 

Torsion -2.7870 + 40.23271i -2.7932 + 39.1 123i -2.7780 + 40.2057i 
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8.3 Validation of Multiple Branch Formulation 

The formulation for multiply-branched blades is validated by 
performing the calculations for the twin beam model considered by 
Sivaneri and Chopra [33]. The data for the twin beam model are 
shown in Table 8.6. 

The rotating natural frequencies are shown in Table 8.7. 
Frequencies about both the initial state and the nonlinearly 
deformed trim state are presented. Those computed by the present 
method agree well with Sivaneri and Chopra's results [33]. 

The nonlinear steady state trim deflections computed from the 
present method are compared with those obtained in [33] in Table 
8.8. Minor differences in the trim deflection are attributable to 
differences in input data in the axial stiffnesses of the two load path 
branches at the blade root. Specifically, the axial stiffnesses of these 
members significantly effects the in-plane bending stiffness due to 
load path offsets. This parameter was not explicitly defined by 
Sivaneri and Chopra in [33] and thus may differ from that used for 
the present calculations. 

The stability eigenvalues for both methods are compared in 
Table 8.9. Examination of the real parts of the flap and torsion 
eigenvalues shows good agreement, with both methods indicating a 
negative sign and subsequently stable flap and torsion motions. 
However both methods result in a lag eigenvalue with a positive real 
part. This indicates that the lag mode of the system is unstable, and 
any perturbation motion will grow with time. The difference in 
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Table 8.6 

Data for the branched blade 



Blade 

2 Identical Load Paths 


C/R = rc/40 

0 = 0.0 


Blade/R = 0.75 

P pc = 0.05 rad = 2.8648° 


Flexure/R = 0.25 

EI y /£2 2 m 0 R 4 = 0.007243 


0 = 12.7183 (constant) 

EI z /f2 2 m oR 4 = 0.083454 


P pc = 0.05 rad = 2.8648° 

GJ/£2 2 m • R = 0.000925 


EI y /Q 2 m 0 R 4 = 0.014486 

km j/R “ 0-0 

— 

EI z /fl 2 m 0 R 4 = 0.166908 

EA/£2 2 m R 4 = 0.0604646 


GJ/Q 2 m 0 R 4 = 0.0004625 

k m2 /R = 0.0125 

— 

m/m 0 R 4 = 1 

k A/km = 0.5 


e/R = 0.0 

e/R = 0.0 


kmr 0.0 
km 2 = 0.025 

m/m Q = 0.5 


(k A /k m ) 2 = 1.5 
EA/£2 2 n 0 R 4 = -0.1209293 

Clevis Geometry for Twinbeam Model 

— 

a = 0.1 

hyjR/C = 16 


a = 6.0 per rad 

h Z j = 0.0 in 

— 

c do = 0.0095 

hy 2 R/C = 16 


0 rc = 0.0° 

h Z2 = 0.0 in 
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Table 8.7 

Natural frequencies (co/£2), twin branched blade 


Initial State Trimmed state 


Ct/cf = 0 Ct/ct = 0.1, p pc = 0.05 (2.86°); y = 5.0 


Present 

Formulation 

Ref. 33 

Present 

Formulation 

Ref. 33 

1.150 

1.149 

1.14888 

1.14974 

1.874 

1.870 

1.77702 

1.77961 

2.908 

2.910 

2.90542 

2.88141 

3.675 


3.78192 


7.624 


7.62226 


8.536 


8.41963 


10.007 


10.14406 


11.233 


11.23863 


15.224 


15.22107 


15.7355 


15.77955 
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Table 8.8 

Nonlinear steady state tip deflections, twin branched blade 
Ct/ct =0.1, Pp C = 0.05 (2.86°); y= 5.0, (o w /Q = 1.15, co v /n = 1.87 
co<t,/£2 = 2.91, zero inboard pitch 


U 0 Wo V 0 


<t>0 


Present TM 0.001548 0.01595 -0.004130 -0.02970 

Formulation 


Ref. 33 


0.01613 


0.011213 


-0.003750 


-0.030735 
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Table 8.9 

Stability eigenvalues, twin branched blade 









Flap (X/Q) 

Lead-lag (X/Q.) 

Torsion (X/Q.) 


Present TM 

-0.34539 + 1.03535i 

0.01102 + 1 .761 09i 

-0.40113 + 2.92296i 


Formulation 





Ref. 33 

-0.35267 + i 

0.00445 + 1.76i 

-0.38119 i 
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magnitude of the real part of the lag eigenvalues can be due to 
differences in the trim position the stability is computed about. 
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9 . CONCLUDING REMARKS 

The objective of the present research is to extend the transfer 
matrix method to a new class of problems and generate numerical 
results to verify the concepts developed. Herein a direct transfer 
matrix method is developed to determine the dynamic characteristics 
of branched autonomous nonlinear rotor blades. The new features of 
the present formulation compared to traditional transfer matrix 
methods are its ability to treat nonlinear boundary value problems 
(using a Newton-Raphson iteration scheme developed for distributed 
systems) and treat multiply branched distributed systems. In the 
case of the multiply branched rotor systems, a rapid iterative scheme 
is employed for the estimation of the tension coefficients in the blade 
root branches. 

The analysis is coded in a FORTRAN computer program, which 
calculates (1) natural frequencies and mode shapes about both the 
initial undeformed and deformed trim states, (2) nonlinear steady 
state deflections corresponding to the hover trim state and (3) 
aeroelastic stability characteristics of single and multiple-branch 
rotor blades. Throughout the actual calculations, the order of the 
matrices involved is only six by six so the method is computationally 
efficient. 

The analysis is applied to two different rotor configurations. A 
conventional single-branch blade is considered to validate the 
nonlinear portion of the analysis. The single-branch blade 
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frequencies, nonlinear trim deflections, and complex stability 
eigenvalues presented provide excellent correlation with the known 
results. A twin-branch blade is also analyzed and its frequencies, 
nonlinear trim deflections and complex stability eigenvalues are in 
agreement with the published data. 

The numerical results thus validate the advancement of the 
transfer matrix method to treat nonlinear distributed boundary 
value problems with multiple branches. 

The extended transfer matrix method has a great potential for 
use in several classes of engineering problems because of its 
computational efficiency. With computer speed roughly doubling 
every eighteen months, it is conceivable to tackle system 
optimization or near real-time system simulation with an 
unprecedented level of modeling sophistication. This could radically 
change the way current designs are developed, because it would 
allow an designer to evaluate many more designs and thus explore a 
much bigger region of the design space. 

The next logical step for the rotorcraft application is to attack 
the forward flight problem. The primary complication in forward 
flight is the changed nature of the governing equations. In 
particular, the aerodynamic lift, drag and moment vary in a periodic 
fashion around the azimuth. This results in a set of nonlinear partial 
differential equations in space and time with periodic coefficients. In 
hover, the equations have constant coefficients and by assuming the 
motions are simple harmonic motions it is possible to reduce the 



partial differential equations in space and time to ordinary 
differential equations in space. 
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This approach is not viable for the forward flight equations. 
There are several alternatives for computing the solution of the 
periodic forward flight equations. 

1 . The free vibration problem is solved and the modes 
(eigenfunctions if using a continuous system model, 
eigenvectors if using a discrete system model) are used to 
reduce the partial differential equations to a set of ordinary 
differential equations in generalized time coordinates. The 
solution is then obtained by integrating the ordinary 
differential equations in generalized coordinates and 
constructing the complete solution via modal 
transformations. In this approach the transfer matrix 
method is only used to solve the free vibration problem. 

2. The second approach is known as harmonic analysis. In this 
approach the rotor blade is not modeled as a continuous 
system but is discretized into a finite set of elements. The 
state vector at a given radial location undergoes periodic 
variation as it moves around the azimuth. Thus the motion at 
that station can be expanded in a Fourier series which has as 
its frequencies multiples of rotor speed. For an infinite 
Fourier series the result is an infinite set of algebraic 
equations, with one set of algebraic equations each for the 
zeroth, first, second, etc. harmonic coefficients in the Fourier 
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series. In practice the Fourier series are truncated and so the 
result is a finite set of algebraic equations which are solved 
with routine linear algebra techniques. This approach lends 
itself readily to the transfer matrix method since transfer 
matrices can be used to relate the harmonic coefficients at 
different radial locations. 

Both approaches are adequate for obtaining solutions to the 
forward flight equations. However when using the modal approach 
some approximations are incurred, and the level of accuracy is 
dependent on how many modes are used. Highly nonlinear systems 
are unwieldy to analyze using modes, and usually require careful 
formulation and a large number of modes for reasonable accuracy. 
The accuracy of the combined Fourier series and transfer matrix 
approach is only dependent on the number of terms retained in the 
series, and nonlinear systems are easier to model. For these reasons 
the second approach is currently under development at Boeing 
Helicopters. It is the foundation for the work in progress described 
in reference [45]. 
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Appendix A: Helpful Integrals to Evaluate the 

Coefficients Defined in 
Equations (3.47), (2.48) 
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